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Sean Phillips and Ricardo G. Sanfelice^ 


Abstract 

The property of desynchronization in an all-to-all network of homogeneous impulse-coupled oscillators 
is studied. Each impulse-coupled oscillator is modeled as a hybrid system with a single timer state that 
self-resets to zero when it reaches a threshold, at which event all other impulse-coupled oscillators adjust 
their timers following a common reset law. In this setting, desynchronization is considered as each 
impulse-coupled oscillator’s timer having equal separation between successive resets. We show that, 
for the considered model, desynchronization is an asymptotically stable property. For this purpose, we 
recast desynchronization as a set stabilization problem and employ Lyapunov stability tools for hybrid 
systems. Furthermore, several perturbations are considered showing that desynchronization is a robust 
property. Perturbations on both the continuous and discrete dynamics are considered. Numerical results 
are presented to illustrate the main contributions. 


1 Introduction 

Impulse-coupled oscillators are multi-agent systems with state variables consisting of timers that evolve 
continuously until a state-dependent event triggers an instantaneous update of their values. Networks of 
such oscillators have been employed to model the dynamics of a wide range of biological and engineering 
systems. In fact, impulse-coupled oscillators have been used to model groups of fireflies [I], spiking neurons 
Elia, muscle cells [4] , wireless networks [5] , and sensor networks • With synchronization being a property 
of particular interest, such complex networks have been found to coordinate the values of their state variables 
by sharing information only at the times the events/impulses occur [TJE]- 

The opposite of synchronization is desynchronization. In simple words, desynchronization in multi¬ 
agent systems is the notion that the agents’ periodic actions are separated “as far apart” as possible in 
time. Desynchronization is similar to clustering or splay-state configurations, and is sometimes referred 
in the literature as inhibited behavior [SlIS]- For impulse-coupled oscillators, desynchronization is given 
as the behavior in which the separation between all of the timers impulses is equal [10) . This behavior 
has been found to be present in communication schemes in fish m and in networks of spiking neurons 
[Hill. Desynchronization of oscillators has recently been shown to be of importance in the understanding 
of Parkinson’s disease [Ml [IS], in the design of algorithms that limit the amount of overlapping data transfer 
and data loss in wireless digital networks [5] , and in the design of round-robin scheduling schemes for sensor 
networks |6]. 

Motivated by the applications mentioned above and the lack of a full understanding of desynchronization 
in multi-agent systems, this paper pertains to the study of the dynamical properties of desynchronization 
in a network of impulse-coupled oscillators with an all-to-all communication graph. The uniqueness of the 
approach emerges from the use of hybrid systems tools, which not only conveniently capture the continuous 
and impulsive behavior in the networks of interest, but also are suitable for analytical study of asymptotic 
stability and robustness to perturbations. 

More precisely, the dynamics of the proposed hybrid system capture the (linear) continuous evolution 
of the states as well their impulsive/discontinuous behavior due to state triggered events. Analysis of the 
asymptotic behavior of the trajectories (or solutions) to these systems is performed using the framework of 
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hybrid systems introduced in [min]. To this end, we recast the study of desynchronization as a set stabi¬ 
lization problem. Unlike synchronization, for which the set of points to stabilize is obvious, the complexity 
of desynchronization requires first to determine such a collection of points, which we refer to as the desyn¬ 
chronization set. We propose an algorithm to compute such set of points. Then, using Lyapunov stability 
theory for hybrid systems, we prove that the desynchronization set is asymptotically stable by defining a 
Lyapunov-like function as the distance between the state and (an inflated version of) the desynchronization 
set. In our context, asymptotic stability of the desynchronization set implies that the distance between 
the state and the desynchronization set converges to zero as the amount of time and the number of jumps 
get large. Using the proposed Lyapunov-like function and invoking an invariance principle, the basin of 
attraction is characterized and shown to be the entire state space minus a set of measure zero, which turns 
out to actually be an exact estimate of the basin of attraction. Furthermore, also exploiting the availability 
of a Lyapunov-like function, we analytically characterize the time for the solutions to reach a neighborhood 
of the desynchronization set. In particular, this characterization provides key insight for the design of al¬ 
gorithms used in applications in which desynchronization is crucial, such as wireless digital networks and 
sensor networks. 

The asymptotic stability property of the desynchronization configuration is shown to be robust to several 
types of perturbations. The perturbations studied here include a generic perturbation in the form of an 
inflation of the dynamics of the proposed hybrid system model of the network of interest and several kinds of 
perturbations on the timer rates. Using the tools presented in [Min], we analytically characterize the effect 
of these perturbations on the already established asymptotic stability property of the desynchronization 
set. In particular, these perturbations capture situations where the agents in the network are heterogeneous 
due to having differing timer rates, threshold values, and update laws. To verify the analytical results, we 
simulate networks of impulse-coupled oscillators under several classes of perturbations. Specifically, we show 
numerical results when perturbations affect the update laws and the timer rates. 

The remainder of this paper is organized as follows. Section [2] is devoted to hybrid modeling of networks 
of impulse-coupled oscillators. Section ITT] introduces an algorithm to determine the desynchronization set. 
Section 13.21 presents the stability results while the time to convergence is characterized in Section 13.31 The 
robustness results are in Section 13.41 Section |4] presents numerical results illustrating our results. Final 
remarks are given in Section [5] 

Notation 

• R denotes the space of real numbers. 

• M" denotes the n-dimensional Euclidean space. 

• N denotes the natural numbers including zero, i.e., N = {0,1, 2,...}. 

• For an interval /C = [0,1] and n £ N \ {0}, K-n is the n-product of the interval /C, i.e., /C„ = [0,1] x 
[0,1] X ... X [0,1]. 

• B is the closed unit ball centered around the origin in Euclidean space. 

• 1 is an N column vector of ones. 

• 1 is an N X iV matrix full of ones. 

• I is the N X N identity matrix. 

• Given a closed set ^ C M" and x £ M”, \x\a '■= min^g^ \x — z\. 

• Given x £ R”, |a;| denotes the Euclidean norm of x. 

• The c-level set of V : domU —>■ R is given by Lv{c) := {x £ domU :V{x) = c}. 
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Figure 1: An example of two impulse-coupled oscillators reaching desynchronization (as Ati converges to a 
constant.) The internal resets ( dark red circles) map the timers to zero. The external resets ( light green 
circles) map the timers to a fraction (1 -|- e) of their current value. 


2 Hybrid System Model of Impulse-Coupled Oscillators 

2.1 Mathematical Model 


In this paper, we consider a model of N impulse-coupled oscillators. Each impulse-coupled oscillator has a 
continuous state (r^ for the i-th oscillator) defining its internal timer. Once the timer of any oscillator reaches 
a threshold (r), it triggers an impulse and is reset to zero. At such an event, all the other impulse-coupled 
oscillators rescale their timer by a factor given by (1 -I- e) times the value of their timer, where e € (—1, 0)0 
Figure [1] shows a trajectory of two impulse-coupled oscillators with states ri and T 2 . In this figure, the dark 
red circles indicate when a timer state has reached the threshold and, thus, resets to zero. The light green 
circles indicate when an oscillator is externally reset and, hence, decreases its timer by (1 + e) times its 
current state. 

According to this outline of the model, the dynamics of the impulse-coupled oscillators involve impulses 
and timer resets, which are treated as true discrete events and instantaneous updates, while the smooth 
evolution of the timers before/after these events define the continuous dynamics. We follow the hybrid 
formalism of [Min], where a hybrid system is given by four objects (C, f,D,G) defining its data: 

• Flow set: a set C C specifying the points where flows are possible (or continuous evolution). 

• Flow map: a single-valued map / : —?> defining the flows. 

• Jump set: a set D C specifying the points where jumps are possible (or discrete evolution). 

• Jump map: a set-valued map G : R^ R^ defining the jumps. 

A hybrid system capturing the dynamics of the impulse-coupled oscillators is denoted as T-Ln '■= {C, f,D,G) 
and can be written in the compact form 


Hat: reR^ 


r f = /(r) T e c 

( r+ G G{t) t G D 


^Cf. the model for synchronization in ^ where e > 0. 


( 1 ) 
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where N S N \ {0,1} is the number of impulse-coupled oscillators. The state of 'Hat is given by 

T := [n T2 ... G Pat := 

The flow and jump sets are defined to constrain the evolution of the timers. The flow set is defined by 


C:=Pn, (2) 

where I := {1,2, and f > 0 is the threshold. During flows, an internal clock gradually increases 

based on the homogeneous rate, oj. Then, the flow map is defined as 

/(r) := ojl Vt € C 

with O’ > 0 defining the natural frequency of each impulse-coupled oscillator. The impulsive events are 
captured by a jump set D and a jump map G. Jumps occur when the state is in the jump set D defined as 

D := (r G Pn : G I s.t. n = f} . (3) 

From such points, the i-th timer is reset to zero and forces a jump of all other timers. Such discrete dynamics 
are captured by the following jump map: for each t € D define G(r) = [giir) 32 (t) ... 9 n{t)]^ , where, 

for each i G /, 

f 0 if Ti = f, Tj. < f Vr G / \ {*} 

5i('r) = ■! {0,Ti(l-be)} if Ti = T 3r G/\ {*} s.t. r,. = f (4) 

[ (1 -b e)Ti if Ti < f 3r G I \ {i} s.t. = f 

with parameters e G (—1,0) and f > 0; for r G D, gi is not empty. When a jump is triggered, the state 

Ti jumps according to the z-th component of the jump map gi. When a state reaches the threshold f, it is 
reset to zero only when all other states are less than that threshold; otherwise, if multiple timers reach the 
threshold simultaneously, the jump map is set valued to indicate that either gi(r) = 0 or g^r) = (1 -b 
is possible. This is to ensure that the jump map satisfies the regularity conditions outlined in Section 
For example, consider the case N = 2 the hybrid system 'Hn = {C, f,D,G) has state given by 


GP 2 := [0,t] X [0,f]. 

The states ri and T 2 are the timers for both of the oscillators. The hybrid system 7^2 has the following data: 

G=P2, 

P2 = ' 






fir) = 

1 

1 

Vr 

gG, 


IP 2 

3z G {1, 2} s.t. 

n = T}, 

G(r) = 

9i 
. 92 

('t) 
('t) _ 

Vr G D, 


are defined as 







if Ti 

= f,T2<f 


' 0 



if T2 = f, 

Tl 

if Ti 

= f,T2=f 

52(t) = < 

{0,r2(l-be 

)} 

if T2 = T, 

Tl 

if n 

<f,T2=f 


, (l+e)r2 


if T2 < f, 

Tl 


= r 
= r 


2.2 Basic Properties of T-Ln 

2.2.1 Hybrid Basic Conditions 

To apply analysis tools for hybrid systems in m, which will be summarized in Section [31 the data of the 
hybrid system Pn must meet certain mild conditions. These conditions, referred to as the hybrid basic 
conditions^ are as follows: 

^In [8], a more general flow map and a jump map incrementing r; by e > 0 are considered. 
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Al) C and D are closed sets in 
A2) / : —>■ R-^ is continuous on C. 

A3) G : R^ =4 R^ is an outer semicontinuoufH set-valued mapping, locally bounded on D, and such that 
G(x) is nonempty for each x G D. 

Lemma 2.1 T-L^ satisfies the hybrid basic conditions. 


Proof Condition (Al) is satisfied since C and D are closed. The function / is constant and therefore 
continuous on C, satisfying (A2). With G as in (jH), the graph of each gi is defined as 

gph(gz) = {[x,y) : y G g^ix),x G D} 

= {{x,y) : y = 0,Xi = f,Xr < f Vr ^ i,x G D} U {{x,y) : y = {1 + e)xi,Xi < f 3xr = r,x G D} 

which is closed. Then the set-valued mapping G is outer semicontinuous. By definition, G is bounded and 
nonempty for each t € D, and hence it satisfies (A3). ■ 

Note that satisfying the hybrid basic conditions implies that 'Hat is well-posed m Theorem 6.30], which 
automatically gives robustness to vanishing state disturbances; see [min]. Section 13.41 considers different 
types of perturbations that 'Hn can withstand. 

2.2.2 Solutions to T-Ln 

Solutions to generic hybrid systems Ti. with state x G R" will be given by hybrid arcs on hybrid time domains 
defined as follows: 

Definition 2.2 (hybrid time domain) A subset S C R>o x N is a compact hybrid time domain if 

j-i 

5'= U (fe,ii-n],j) 

j=o 

for some finite sequence of times Q = t^ < ti < t^ ■■■ < tj. A subset S C R>o x N is a hybrid time domain 
if for all (T, J) G S, S H ([0,T] x {0,1, ■■■J}) is a compact hybrid time domain. 


Definition 2.3 (hybrid arc) A function x : doma; —>■ R" is a hybrid arc if domx is a hybrid time domain 
and if for each j G N, the function 1 1—^ x(t,j) is locally absolutely continuous. 

Definition 2.4 (solution) A hybrid arc x is a solution to the hybrid system H if x{0,0) € CU D and: 

(51) For all j G N and almost all t such that (t,j) G domx, 

x{t,j) G C, x{t,j) = f{x{tj)) . 

(52) For all {t,j) G domx such that (t,j -I- 1) G domx, 

x{t,j) G D, x{tj -P 1) G G{x{t,j)) . 

set-valued mapping G : ^ is outer semicontinuous if its graph {{x,y) : x € € G(a:)} is closed, see |16l 

Lemma 5.10] and m- 
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A solution X is said to be nontrivial if doma; contains at least one point different from (0,0), maximal if 
there does not exist a solution x' such that a; is a truncation of x' to some proper subset of doma;', complete 
if doma; is unbounded, and Zeno if it is complete but the projection of doma; onto ]R>o is bounded. 

Lemma 2.5 From every point in CVJD, there exists a solution and every maximal solution to "Hat is complete 
and bounded. 


Proof The result follows from Proposition 2.10 in [TB] using the following properties. For each point such 
that r € (7, the components of the flow map / are positive and induce solutions that flow towards D. For 
each T G D, the jump map satisfies G{t) C C. Since it is impossible for solutions with initial conditions 

r(0,0) € C U H to escape CUD, all maximal solutions are complete and bounded. ■ 

Due to the jump map G, if the elements of the solution are initially equal (denote this set as S := {t € Pn '■ 
Bi,r € I,i r,Ti = Tr}) it is possible for them to remain equal for all time. Furthermore, it is also possible 
for solutions to be initialized on the jump set such that one element is at the threshold and another is equal 
to zero then after the jump they will be equal, e.g. let ti = f, T2 = 0 then = 0. We denote this set 

as Q := {t £ D \ S : 3i, r G I, i r, Ti = 0, Tr = f}. The next result considers solutions initialized on the set 

X :=sug. 

Lemma 2.6 For each t(0, 0) € X, there exists a solution r to T-Ln from r(0,0) such that, for some M £ 

{0, 1}, T(t,j) £ S for all t + j > M, (t,j) £ domr. 


Proof Consider a solution r to the hybrid system Hn with initial condition r(0, 0) £ S. Due to the flow 
map for each state being equal, r remains in S during flows. Furthermore, at points r £ SUD, the jump map 
G is set valued by the definition of in (jH). From these points, G(t) fl 5 ^ 0. In fact, for each r(0, 0) £ S, 
there exists at least one solution such that T{t,j) £ S for all t + j > 0, with {t,j) £ domr. Consider the case 
of solutions initialized at t(0, 0) £ Q (Note that r(0,0) £ D). It follows that for some r £ I, Tr{0, 0) = f and 
gr{T{0,0)) = 0. Therefore, after the initial jump, we have that G(t( 0,0)) C 5 ^ 0, by which using previous 
arguments implies that T{t,j) £ S for all t + j > I. 

Furthermore, there is a distinct ordering to the jumps. If r is such that Ti ^ Tr for a\\ i ^ r then the 
ordering of each Ti is preserved after N jumps. More specifically, we have the following result. 

Lemma 2.7 For every solution t to "Hat with r(0,0) X, if at (tj,j) £ domr we have 

0 < < Ti^{tj,j) < ... < Ti^{tj,j) < f 

for some sequence of nonrepeated elements {im}m=i of I (that is, a reordering of the elements of the set 
I = {1,2,..., N}) then, after N jumps, it follows that 


0 < T,^{tj+Nj + N) < n^{tj^N,j + N) < ... < Ti^{tj^N,j + N) <T. 


Proof Let r be a solution to Hn from Pa; \ A". There exists a sequence ik of distinct elements with ik £ I 
for each k £ I, such that 0 < Ti„{t,j) < Ti,^{t,j) < ... < Tij.,(t,j) < f over [tojG] x {0}. After the jump 
at {t,j) = (ti,0) we have 0 = Ti^(t,j + 1) < n^{t,j + I) < Ti.,{t,j + 1) < ... < n^_,{t,j + 1) < f. 
Continuing this way for each jump, it follows that after iV — I more jumps, the solution is such that 
0 < Ti^{tN,j + N) < Ti 2 {t]Y,j + N) < ... < Tij^{tM,j + N) < f and the order at time {t,j) is preserved. 

Using these properties of solutions to Hn, the next section defines the set to which these solutions converge 
and establishes its stability properties. 
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3 Dynamical Properties of 7/jv 


Our goal is to show that the desynchronization configuration of T-Ln, which is defined in Section |3.11 is 
asymptotically stable. We recall from [min] the following definition of asymptotic stability for general 
hybrid systems with state x G M". 

Definition 3.1 (stability) A closed set A C K" is said to be 

• stable if for each e > 0 there exists 5 > 0 such that each solution x with |a:(0,0)|^ < 5 satisfies 
\^ifA)\A ^ ^ /c”' (fA) S domx; 

• attractive if there exists /r > 0 such that every maximal solution x with |a;(0,0)|^ < fi is complete and 
satisfies 

Ihnp ^ 0 ? 

• asymptotically stable if stable and attractive; 

• weakly globally asymptotically stable if A is stable and if, for every initial condition, there exists a 

maximal solution that is complete and satisfies limp j)g(jomx,t+j^oo = 0 . 


The set of points from where the attractivity property holds is the basin of attraction and excludes all 
points where the system trajectories may never converge to A. In fact, it will be established in Section [321 
that the basin of attraction for asymptotic stability of desynchronization oi Hn does not include any point 
T such that any two or more timers are equal or become equal after a jump, which is the set A defined in 
Lemma 1221 For example, consider the case N = 2, i.e., 7 ^ 2 - Then, the set X 2 is defined as 


^”2 = ^2 U 172 = ({r £^ 2 :^ 1 = T2}) U{{t G D : gi{T) = T2} U {t G D : g 2 {T) = ri}). 


(5) 


Note that the set S 2 defines the line n = T 2 in P 2 and G 2 is given by the points {(0, r), (r,0)} in P 2 ; see 
Figure 2(a) For N = 3, the set X 3 is defined as 


<^3 — <$3 U G 3 


( 6 ) 


where 


^3 ={t G P3 : Ti = T2} U {r G P3 : n = T3} U {r G P3 : T2 = T3} 


(7) 


and 


G3={tGP3 ■■ gi{T) = T3} U {t G P3 : ff 2 (T) = Ti} U {t G P3 : 53(t) = n} U {r G P3 :ff2(T)=T3} 
U{tGP3 : 53(t) = T 2 } U {t G P3 :gi(T)=T 2 }. 

Then, T 3 is defined by the union of <$ 3 , which is the as the union of the planes in P 3 given by ri = T 2 , 
Ti = T 3 , and T 2 = T 3 , and G 3 , which is given by {(r, r 2 , 0 ), ( 0 , r 2 , f), (n, f, 0 ), 

(ti, 0, r), (f, 0, r 3 ), (0, f, Ts) : r G P 3 }; see Figure [3] For this purpose, a Lyapunov-like function be con¬ 
structed in Section 13.21 to show that a compact set denoted A, defining the desynchronization condition, is 
asymptotically stable and weakly globally asymptotically stable. 

3.1 Construction of the set A for T-Lpf 

In this section, we identify the set of points corresponding to the impulse-coupled oscillators being desynchro¬ 
nized, namely, we define the desynchronization set. We define desynchronization as the behavior in which 
the separation between all of the timers’ impulses is equal (and nonzero), see Figured] More specifically 
desynchronization is defined as follows: 



Tl 




(a) The sets X 2 (dashed green) and A (solid blue and (b) A simulation (dashed blue) of "^2 showing the at- 
red). The set A is defined by the union of and I 2 tractivity of the set A (solid black). 

(See Section l3.1|l . 


Figure 2: Sets associated with 'H 2 and a solution to it from t(0,0) 


[0.7,0.75]^ with e: = 0.3 and r 



Figure 3: Points in each blue plane and line belong to the set X^. 
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Definition 3.2 A solution t to 'Hn is desynchronized if there exists A > 0 and a sequence of non-repeated 
elements {*m}m=i 0 / ^ (that is, a reordering of the elements of the set I = {1,2,..., iV}^ such that 
limj_>oo(f}”* - = A for all m G {1,2,..., N - 1} and limj^,^{t^ - tj^) = A, where 

the sequence of jump times of the state Ti^. 


In fact, this separation between impulses leads to an ordered sequence of impulse times with equal separa¬ 
tion. The desynchronization set A for the hybrid system T-i^ captures such a behavior and is parameterized 
by £, the threshold f, and the number of impulse-coupled oscillators N. 

To define this set, first we provide some basic intuition about the dynamics of Ttw when desynchronized. 
The set A must be forward invariant and such that trajectories staying in it satisfy the property in Defini¬ 
tion o Due to the definition of the flow map /, there exist sets in the form of “lines” ik, each of them 
in the direction 1, which is the direction of the flow map, intersecting the jump set at a point which, for 
the fc-th line, we denote as r^. We define the desynchronization set as the union of sets collecting points 
T = + \s G Pn parameterized by s S M. Figure | 2 (a)| shows and ^2 (solid blue and red) for the case 
N = 2. 


To identify r^, consider a point G D\X with components satisfying rf = f > > ... > t%. 

Due to Definition [321 it must be true that the difference between jump times are constant. This means that 
there must be some correlation between A and the difference between, in this case, rf and Moreover, 
there must be a correlation between rf and all other states at jumps. It follows that this point belongs to 
A only if the distance between the expiring timer (rf) and each of its other components (rf', f€/\{l})is 
equal to the distance between the value after the jump of the timer expiring next (tJ and the value after 
the jump of its other components (rf ■*", i G I \ {2}), respectively. This property ensures that, when in the 
desynchronization set, the relative distance between the leading timer and each of the other timers is equal, 
before and after jumps. More precisely, 


T, - T,- = T. 




-k + 


xt(i) 


V*e J\{1}, 


(9) 


where = G(t^) and next(i) = i-|-lifi-|-l<A and 1 otherwise^ Since X contains all points such 
that at least two or more timers are the same, we can consider the case when one component of is equal 
to f at a time. For each such case, we have {N — 1)! possible permutations of the other components and N 
possible timer components equal to f, leading to N\ total possible sets ik- 

To illustrate computation of in and the construction of A, consider the case of A = 2 and 
t} = f > To. For f = 2, (O becomes 

^-^2 =^2(£ + 1) 

which leads to rf = It follows that = [f, Similarly for rf = f > t(, we get from (|9]) 

the equation f — t( = rfle -I- 1), which implies ■ A glimpse at the case for A = 3 with 

tI = f >Tf > T 3 indicates that m leads to 


t-tI= tI(1 -be) - T 3 i(l -be), t - T 3 ^ = rj(l -b e) - 0. 


The solution to these equations is A = [t, r(e -b 2)/(e^ -b 3e -b 3), r/(e^ -b 3e -b 3)]'’". 

For the A case, the algorithm above results in the system of equations Fr^ =h, where 


F = 


0 


0 


0 ( 2 -be) -( 1 -be) 

0 (I-be) I 

0 (I-be) 0 

0 (I-be) 0 


0 

0 

-(l+£) 

1 

0 

0 


0 


-(1-be) 

1 


(10) 


^Note that G is single valued at each A X. 
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and b = fl, where Tg is the state sorted into decreasing order. For example, if is such that = f > 
ff > Tg , then Ts is given as [tJ, , Tg ]^. It can be shown that for any e S (—1,0), a solution Tg exists (see 
Lemma Fa. II) . Then, Tg needs to be unsorted and becomes in the definition of the set £k- 

The solution to Tts = b is the result of a single case of t S D \ A. As indicated above, to get a full 
definition of the set A, the A^! sets £k should be computed. For arbitrary TV, the set A is given as a collection 
of sets £k given by 

N\ 

fe=i 

where, for each fc S {1, 2,..., A^!}, := {t : t = + Is G Pat, s £ K}. For the case N = 2, the points 

for k G {1,2} lead to the set A given by 


A — £iU £2 



+ Is G P2, s G 


U 



r = 


e+2 

r 


+ Is G P2, s 



Figure 2(a) shows these sets in the (rg,T 2 )-plane (solid blue and red). Figure 2(b) shows a solution to 'H 2 - 
The initial conditions for the simulation are t(0, 0) = (0.75, 0.7). 

Furthermore, for the case N = 3 the points for k G (1, 2,..., 6 } lead to the set ^3 given by 


A3 = £1 U £2 U £3 U £4 U £5 U ^6 



r 




r 


II 

(e+2)T 

e2+3e+3 

T 

H- Is ^ -fsjS ^ M 

M 

II 

r 

e^+3e+3 

(e+2)i- 

-|- Is G -P3, s G M 1 


L i^+3e+3 J 

j 



L e^-|-3e+3 J 

) 


r : T = 

r (£+ 2 ) 1 - 1 

T 

+ Is € -P 3 , s € M. 

hi 

T : T = 

r T -| 

e^+3e+3 

f 

H- Is G P 3 -)S G M. , 

1 

: r = 

T 

L e^-|-3e+3 J 
r {e+2)f n 
e^^T3e+3 

T 

J 

-f Is G P 3 , s G M 

I 1 

K-r- 

(e+2)t 

L e2+3e+3 J 

r T -| 

e24-3£:+3 

(£+2)t 

i 

-|- Is G P 3 ,s G 

e^+3e+3 

r 


e2+3e:+3 

f 


Figure 4(a) shows these sets in the (n, T 2 , T 3 )-plane (solid colored). Figure [4(b) | shows two solutions to 
T-Ls- Note how each simulation has jumps that take the trajectory close to different lines. This is due to 
a preservation of order for each Ti as seen in Lemma 12.71 This preservation of order will be used in the 
Lyapunov stability proof in the next section. 


3.2 Lyapunov Stability 

Lyapunov theory for hybrid systems is employed to show that the set of points A is asymptotically stable. 
Our candidate Lyapunov-like function, which is defined below and uses the distance function, is built by 
observing that there exist points where the distance to A may increase during flows. This is due to the sets 
£k being a subset Pn. To avoid this issue, we define 

N\ 

A=[j£k^A 

k^l 

where ik is the extension of ik given by 

4 = {r gK^ :T = T'=-bls,sGK}. ( 12 ) 
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(a) The desynchronization set A for N = 3, with 
f = 1 and e = —0.3. 



(b) Two solutions (magenta and cyan) to l-Ls such 
that r(0, 0) ^ with e = —0.3 and f = 1 showing 
the attractivity to As = 


Figure 4: (a) Set of points {h, (.2 ,■■■,(■&) defining As and |(b)| two simulations (dashed, one in cyan and the 
other in magenta) converging to the set ^3 (solid colored). 


Then, with this extended version of A, the proposed candidate Lyapunov-like function for asymptotic stability 
of .4. for 'Hat is given by the locally Lipschitz function 

1/(t) =min{|r|j^,|T|j^,...,|r|j^,...,|r|j^,} V r e Pv \ 4” (13) 

where, for some k, \t\j^ is the distance between the point r and the set ffcH The following theorem establishes 
asymptotic stability of ^ for Hat. We show that the change in V during flows is zero and that at jumps we 
have a strict decrease of V; namely, V(G(t)) — V(t) = —\e\V{T). A key step in the proof is in using [HI 
Theorem 8.2] on a restricted version of T-Ln- 

Theorem 3.3 For every TV £ N, TV > 1, f > 0,a; > 0, and e £ (—1,0), the hybrid system T-Ln is such that 
the compact set A is 

1 . asymptotically stable with basin of attraction given by B_a, := Pjq \ X. 

2 . weakly globally asymptotically stable. 


Proof Let the set Xy define the u-inflation of X (defined in Lemma 1231) . that is, the open sell^ Xy := {r £ 
: \t\x < where v £ (0,u*) and v* = min^g^,^, ja; — y\. Given any v £ (0,u*), we now consider a 

restricted hybrid system Hn = {f,C,G,D), where G := G \ Xy and D ■.= D \ Xy, which are closed. We 
establish that A is an asymptotically stable set for T-Ln- _ 

Note that the continuous function V, given by m, is defined as the minimum distance from t to A, 
where A is the union of A^! sets ik in ([T^ . To determine the change of V during flowl3, we consider the 
relationship between the flow map and the sets ik- The inner product between a vector pointing in the 
direction of the set ik and the flow map on G satisfies 

=ojN = |l||wl| = |1 ||/(t)|cos6» 

®The set can be described as a straight line in M"' passing through a point and with slope 1. Then, can be written 

as the general point-to-line distance |(r^ — r) — 1/A^((r^ — 

®The set Xy is open since every point r E Xy is an interior point of Xy. 

^ Its derivative can be computed using Clarke’s generalized gradient m- 
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, which is only true if 6 is zero. Therefore, the^direction of the flow map and of the vector defining are 
parallel, implying that the distance to the set A is constant during flow^ ^ 

The change in V during jumps is given by V{G(t)) —V(t) for r G D \ A. Due to the fact that we can 
rearrange the components of r G \ without loss of generality, we consider a single jump condition, 
namely, we consider r such that f = ti > T 2 > ... > tn-i > Using the formulation in Section Id.ll and 

Lemma [01 the elements of the vector associated with for this case of r are given by rf = ^W-i) - t-t, 

which by Lemma fA.2l is equal to After the jump, G(t) is single valued and is such that its 

elements are ordered as follows: 52 ('r) > gsir) > ... > gN{T) > 51 (t) = 0. Specifically, the jump map is 
G(t) = [0, (1 + e)T 2 ,..., (1 + ejrjv]^. Then, the formulation in Section ITT] and Lemma [AT] leads to a case 
of denoted as . By Lemma IA.21 the elements of the vector are given by r^' = and 

rf = for i > 1. Due to the ordering of r and G(t), is a one-element shifted (to the right) 

version of . 


From the definition of above, U at r reduces to 


= klj, = 




for some k. Note that 


N 


N 


(r^--r)^l = ^T'=-^r, 
i—1 i—1 

reduces to ~ Y^f =2 P since ti =t^ = f. Using Lemmas lA.21 and lA.31 it follows that 

Et2Epl-;(£ + iF__ ((£ +1)^-1)-iVe_ 




EEo'(^ + 1 )" 


-r = 


e((e -I-1)^ - 1) 


-T. 


i=2 

Then, the first element of the vector inside the norm in the expression of V (r) is given as 

^ 7 V N ' 




1 / ((e-bl)^-l)-A^e_ 


e((e -I-1)^ - 1) 


-E' 

2=2 


N 


{{e + l)^ - 1) - Ne 1 
£A^((e +1)^-1) 


2=2 


while the elements with are given by 


- y 


2 = 2 


e((e +1)^-1) 


(£-b 1 )^-’"+! - 1 . 
{e+1)^-1 


-T -Trr, \- 


1 / ((£ -b 1)^ - 1) - iVe 


N 


N 


£((£ -b 1)^ - 1) 




eN{e + - ((e + l)^ - l) 

-r 


£iV((£-b 1)^ - 1) 


N-1 

N~' 


N 




2 = 2 , 2^771 


After the jump at r, since G(r) is single valued, U(G(r)) is given by 


|G(r)|,y = 


(7'='-G(r))-y((?^-'-G(r))^l)l 


1 


-k' 
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Note that {t'^' - G(t))^ 1 = 5i('^) reduces to “ ^^ 2 ( 1 +^)'^*’ since gi(r) = 0 and 

giir) = (1 + e)Ti for i > 1 . Using Lemmas lA.21 and lA.3[ it follows that 


N 



i=l 


El^E^Joie+iy _ 

(£ + l)((£ + l)^-l)-jV£_ 

£((£ + 1 )^ - 1 ) 


which leads to 


(f‘' - G{r)fi = 


£((£ + 1 )^ - 1 ) 

The first element inside the norm in V{G{t)) is given by 


(rf -3i(r))-^ 


1 / (£ + 1 )((£ + 1 )" — 1 ) — Ne 


£((£ + 1 )^- 1 ) 


N 




£ _ (£ + 1)((£ + 1)^ — 1) —-/V£ _ 

-T - 7 - 77 ——rnrr - 77 - t 


(£ + 1)^-1 £iV((£ + 1)^ - 1) 


N 


-J2il+e)n 


V ((£+1)^-1)-N£ 1 ^ 


For each element m > 1, it follows that 


{T^-9m{T)) - 


1 / (£+l)((£ + l)^ -l)-iV£ 


N 


N 


e{{e + 1 )^ - 1 ) 




i =2 


(£ + l)N-m+2 _ I _ 


(^+ 1 ) 


N 


— —r-(l + £) 


fV-1 (£ + l)((£ + l)^-l)-fV£_ 1 


N 


N 


~'Tm 


eN{{e + 1)^ - 1) 


T + 


N 


^ (1 + e)Ti 


N 


£iV(£+l)^-™+l- ((£ + 1)^-1)^ iV_l_ ^ 

- (1 + £)| -nUTTUT-TUv-^- ^ Iv 


£fV((£ + 1)^ - 1) 


N 




Combining the expressions for each of the elements inside the norm of V (G(r)), it follows that V (G(r)) = 
(l + £)U(r). 

Then, the change during jumps is given by U(G(t)) — V(t) = £U(t) where £ G (—1,0). With the 
property of V during flows established above, the change of V along solutions is bounded during flows 
and jumps by the nonpositive functions and ug, respectively, defined as follows: Uq{z) = 0 for each 
z G C and u^{z) = —00 otherwise; ujj[z) = £U( 2 :) for each z G D and uj^{z) = —00 otherwise. Using 
Lemma O the fact that G and D are closed, and the fact th^ every maximal solution to % is bounded 
and complete, by m Theorem 8.2], every maximal solution to approaches the largest weakly invariant 
subset of Lv{r') fl G Cl U (Lug(O) fl G(L„g(0)))] = Lv{r') fl G for r' G V{C). Since every maximal 

solution jumps an infinite number of times, the largest invariant set is given for r' = 0 due to the fact that 
U(G(t)) — U(t) = £U(t) < 0 if r' > 0. Then, the largest invariant set is given by Ly(0) fl G = ^ n G 
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which is identically equal to A. Hence, the set A is attractive. Stability is guaranteed from the fact that V 
is nonincreasing duringjlows and strictly decreasing during jumps. Then, the set A is asymptotically stable 
for the hybrid system "Hjv- We have that A is (strongly) forward invariant and from Theorem 3.4 we know 
that A is uniformly attractive from a neighborhood of itself. Then by Proposition 7.5 in m. it follows that 
A is asymptotically stable. ^ 

Note that the set of solutions to TLn coincides with the set of solutions to TLn from Pn \ Ay. Therefore, 
the set A is asymptotically stable for "Hat with basin of attraction = Pn \ Since v is arbitrary, it 
follows that the basin of attraction is equal to Pn \ X. 

Note that the jump map G, at points r € Tf, is set valued by definition of gi in (j4)). From these points 
there exist solutions to Pn that jump out of X. In fact, consider the case t € X. We have that Ti = Tr for 
some i,r G I. Then, after the jump it follows that 5 i(r) £ {0, (1 + e)f} and grir) £ {0, (1 + ejr}, and there 
exist gi and gr such that gi = gr or gi ^ gr- Since for every point in X there exists a solution that converges 
to A and also a solution that stays in X, X is weakly forward invariant!! 

3.3 Characterization of Time of Convergence 

In this section, we characterize the time to converge to a neighborhood of A. The proposed (upper bound) 
of the time to converge depends on the initial distance to the set A and the parameters of the hybrid system 
(e:,T). 

Theorem 3.4 For every N £ N, N > 1, and every Ci,C 2 such that c > C 2 > ci > 0 with c = 
every maximal solution to Pn with initial condition r(0, 0) £ {Pn \ X) n Ly{c 2 ) is such that 

T{t,j) £ Lv{ci) V(t,j) G domTp + j > M, 

where 



and Lvig) '■= {t £ C VJ D ■. V{t) < g\. 

Proof Let tq = t( 0,0) and pick a maximal solution r to Pn from tq. At every jump time {tj,j) € domr, 
define gi = r(ti, 1), 52 = r(t 2 ,2 ),... ,gj = T{tj, J), for some J G N. From Theorem 13.31 we have that there 
is no change in the Lyapunov function during flows. Furthermore, we have that for each t £ D \ A the 
difference V{G{t)) — V{t) = £V{t) with e G (—1, 0). Since, for every j, T{tj,j) G D, we have 

Viai) -Viro) = eH(ro), 

which implies 

^(51) = (1 + £) V { tq ). 

At the next jump, we have 

V{g2) = (l + e)H(5i) = (1 +e)V(ro). 

Proceeding in this way, after J jumps we have 

y{9j) = (l + e)l^(5J-i) = (1+e)‘^P(ro). 

®For example, consider the case N = 2. If r(0, 0) = [T,f]^ G then there are nonunique solutions due to the jump map 
begin set valued. It follows that after the jump, each ti can be mapped to any point in {0, Ti(l + e)}, which leads to any of 
the following four options of the states (ri, T 2 ) after such a jump: (0, 0), (0, f(l + £)), (f(l + £), 0) or (f(l + e), f(l + e)). If the 
state is mapped to either (0, 0) or (f(l + £), •f(l + e)), then it remains in X 2 . Conversely, if any of the other options are chosen, 
then (ri, T 2 ) leaves X 2 and converges to v4 asymptotically. 
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Figure 5: Time to converge (over r + 1) as a function of £ G [—0.9,—0.1], with C 2 = 0.99r and ci G 
{0.5f,0.3f,0.1f,0.05f} 


From V{gj) = (1 + £)‘^y(ro), we want to find J so that V{gj) < ci when V[tq) < C 2 - Considering the worst 


cast for we want (1 + e)'^C 2 < ci, which implies ^ i therefore J = 


log If 
logTT? 


> 0. For 


each J, the time between jumps satisfies ti — to < —,t2 — ti < . ,tj — tj-i < — . Then, we have that after 


J jumps, J2i=i ^3 ~ ^i-i — Jzi- With to = 0, the expression reduces to tj < = 


log If 
logTT? 


If. Then, after 


t + j > tj + J, the solution is at least ci close to the set A. Defining M = tj + J we then have 


M 



logff 


log 


l+£ 


Figure [5] shows the time to converge (divided by ^ + 1) versus e with constant C 2 = 0.99f and varying 
values of ci. As the figure indicates, the time to converge decreases as |£| increases, which confirms the 
intuition that the larger the jump the faster oscillators desynchronize. 


3.4 Robustness Analysis 

Lemma o establishes that the hybrid model of N impulse-coupled oscillators satisfies the hybrid basic 
conditions. In light of this property, the asymptotic stability property of A for T-Ln is preserved under 
certain perturbations; i.e., asymptotic stability is robust [16]. In the next sections, we consider a perturbed 
version of Hat and present robust stability results. In particular, we consider generic perturbations to 'Hat, 
and two different cases of perturbations only on the timer rates to allow for heterogeneous timers. 

3.4.1 Robustness to Generic Perturbations 

We start by revisiting the definition of perturbed hybrid systems in |16j . 

Definition 3.5 (perturbed hybrid system |16L Definition 6.27]) Given a hybrid system T-L and a func¬ 
tion p : > K> 0 j the p-perturbation of T-L, denoted T-Lp, is the hybrid system 

f X G Cp X G Fp(x) 

\ X G Dp x'^ G Gpix) 
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where 


Cp = {a; e R" : (x + p(x)B) n C 7^ 0}, 

Fp{x) = coni^((x + p(x)B) 0(7)+ p(x)B Vx S R", 

Dp = {x G R" : (x + /o(x)B) n D 7 ^ 0}, 

Gp(x) = {x G R” : X e g + p{g)n, g G G((x + p(x)B) n D)} Vx G R”. 

Using this definition, we can deduce a generic perturbed hybrid system modeling N impulse-coupled oscilla¬ 
tors. Then, for the hybrid system "Wjv, we denote 'Hn,p as the p-perturbation of T-Ln- Given the perturbation 
function p : R^ —)• R>o, the perturbed flow map is given by 

Fp(r) = wl -I- p(r)B ^ t G Cp, 

where the perturbed flow set Cp is given by 

Gp = {r G R^ : (r + p(t)B) n P^v ^ 0}- 

For example, if = 2 and p{t) = p > 0 for all t G R^, which would correspond to constant perturbations 
on the lower value and threshold, then Gp = G -I- pB. The perturbed jump map and jump set are defined as 

Pp = {r G R^ : (r + p(t)B) n P 7 ^ 0}, 

Gp = [pi,p(r),...,p7v,p(r)]^, 

where gi^p is the i-th component of Gp. The following result establishes that the hybrid system T-Lj^ is robust 
to small perturbations. 

Theorem 3.6 (robustness of asymptotic stability) If p : R'^ —^ R>o is continuous and positive on R'^ \ A, 
then A is semiglobally practically robustly ICC asymptotically stable with basin of attraction = P^r \ X, 
i.e., for every compact set K C and every a > 0, there exists S G (0,1) such that every maximal solution 
T to TLn.Sp from K satisfies |T(t, j)|^ < /3 (|t(0, 0)|^, t + j) + a for all (t, j) G domr. 


Proof From Lemma 12.11 the hybrid system Hn satisfies the hybrid basic conditions. Therefore, by [T51 
Theorem 6.8] Hn is nominally well-posed and, moreover, by [161 Proposition 6.28] is well-posed. From 
the proof of Theorem 13.31 we know that the set A is an asymptotically stable compact set for the hybrid 
system Bn with basin of attraction P_4. Since by Lemma l2.51 every maximal solution is complete, then [161 
Theorem 7.20] implies that A is semiglobally practically robustly ICC asymptotically stable. 

Section (4.2.11 showcases several simulations oIBn with p-perturbations on the jump map. 

3.4.2 Robustness to Heterogeneous Timer Rates 

We consider the case when the continuous dynamic rates are perturbed in the form of 

= c(t,j) 

for a given solution r. For example, consider the perturbation of the flow map given by 

/(r) = wl -I- Aw (14) 

where Aw G R" is a constant defining a perturbation from the natural frequencies of the impulse-coupled 
oscillators. Then for some k, during flows, along a solution r such that over [tj,tj+i]><{j} satisfies V{T{t,j)) = 
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\T{t,j)\j^, it follows that c reduces to c{t,j) = Furthermore, the norm of the hybrid 

arc c can be bounded by a constant c given by 


c = 




(15) 


Building from this example, the following result provides properties of the distance to A from solutions r to 
'Hn under generic perturbations on / (not necessarily as in (flS l. 


Theorem 3.7 Suppose that the perturbation on the flow map of T-Cn is such that a perturbed solution r 
satisfies, for each j such that {t : {t,j) G domr} has more than one point, = c{t,j) for all 

t G {t : {t,j) G domr} and T(t,j) G Pn \<T for all {t,j) G domr, for some hybrid arc c with dome = domr. 
Then, the following hold: 

• The asymptotic value o/|r(t,j)|^ satisfies 

lim |r(t,})|j< lim ^(l+e)^"*/ c{t,j)dt (16) 

t-\-J^OO t+J—>-CX) ^ 


• If there exists c > 0 such that |c(t,j)| < c for each {t,j) G domr then 


lim 

t+j—foo 




CT 

klw’ 


(17) 


• If j '■ K>o —^ N *s a function that chooses the appropriate minimum j such that (t,j) G domr for each 
time t and 1 1 —^ c{t,j{t)) is absolutely integrable, i.e., 3B such that 



\c{t,j{t))\dt < B, 


(18) 


then 


lim 

t+j—fco 




B 

e ’ 


(19) 


Proof Consider a maximal solution t to Bn with initial condition r(0,0) G Pn \ X- This proof uses the 
function V from the proof of Theorem 13.31 With V equal to the distance from r to the set A, then, for each 
T G D\X, we have that V{G{t)) — V{t) = sV{t). Using the fact that U(r) = |r|_^ and the fact that, G 
along the solution is single valued, it follows that \t\^ after a jump can be equivalently written as 

+ l)|j = (1 + e)|r(tj, j)|_^. 

By assumption, in between jumps, the distance to the set A is such that ^l'r(C j)l^ = c{t,j), which implies 
that at tj+i the distance to the desynchronization set is given by 


\T{tj+lG)\A = 


c{s,j)ds + \T{tj,j)\j. 


® Let ri^ir) be the vector defined by the minimum distance from r to the line Then, it follows that 

V{t) = ('?■)) ^ • To determine its change during flows, note that on C* \ (A' U A) the gradient is given 

by VV(t) = (rJjT)re^{T)Y = - is given by = 

Sf (ft*" “ ~ W Edilri" - Ti)^) = [X,X,...,X,-1 + X,X,...,X] - the term -1 + ^ corresponds to the j-th el- 

ement of the vector. It follows that (t) = -^1 — I. Then, for each r E C \ (W(r), /(r)) = ( - j /(r). 
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It follows that 


k(ii:0)U= / c(s,0)ds + |t(0,0)|^ 

k(^i, 1)U = (1 + ^) (^J c(s, 0)ds + |t( 0, 0)1 = (1 + e) J c(s, 0)ds + (1 + e)|T(0, 0)| j 

k(i 2 ,l)|^= / c(s,l)ds +(1 + e) / c(s,0)ds +(l + e)|r(0,0)|^ 

Jti Jo 

k(^ 2 , 2)|_4 = (1+ £) ^^ c{s,l)ds + {l+e) J c(s,0)ds + (l + e)|r(0,0)|_^^ . 

Then, proceeding in this way, we obtain 

= (1 + e)^k(o>o )|_4 

. . /■*«+! 

+ ^(1+ £)■’“*/ c{s,i)ds. 

i =0 

For the case of generic t^+i > t > tj, we have that 

(l + £)^k(0,0)|j + ^(l+£)^"* f c{s,i)ds. 

1=0 

Since, we know that as either t or j goes to infinity, j or t go to infinity as well, respectively. The expression 
reduces to 

j pt 3 pt 

lim \T{t,j)\r= lim (1 + £)'^|r(0,0)1 T + lim y^fl + s)-^"* / c(s,i)ds = lim + / c(s,i)ds. 

t-\-j—>-00 J—>-00 J± t-\-j—¥00 1+ 

( 20 ) 


If c{t,j) < c, it follows that 


lim |T(t,j)|^= lim ' f c(s,*)d 


2=0 

3 


< lim y^(l + £)' 

t-\-j—¥oo • ^ 

2=0 


I cdt 


< — lim (1 + ey ’ 
2=0 

cf (1 + ey - 1 

= — hm —--- 

W t+j-s-oo (1 + £) — 1 

cf 

~ \e\uj' 


Lastly, since this hybrid system has the property that for any maximal solution r with (t,j) G domr, if 
t approaches oo then the parameter j also approaches oo, the expression given by limt+j_>oo \T(t,j)\j[ can 

be simplified. To do this, we know that the series + £)'^“* = —— approaches as j —>■ oo. 

Since 1 + £ > 0 for £ S (—1, 0), the series is absolutely convergent and its partial sum sj = 

is such that is a nondecreasing sequence (for each m). This implies that Sj < 1/|£| for all j and for 
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(a) Solutions to 1-12 with r(0, 0) G X 


(b) Solutions to Ms with r(0, 0) G Az 


Figure 6: Solutions to Hat with N G {2, 3} that are initially in the set A. 


each m. Then, it follows that (1 + ey ® for every j,i G N. Since the expression is a function of j only 
and, for complete solutions, t is such that as f —> oo, then j —^ oo, we obtain 


Inn y^(l+£) 

t+i^oo 

i=0 





c(s,i)ds 
|c(s, i)\ds 


< 



|c(s, i)|(is 




c{sA{s))\ds. 


4 Numerical Analysis 

This section presents numerical results obtained from simulating 'Hn- First, we present results for the 
nominal case of given by ©• Then, we present results for Hn under different types of perturbations. 
The Hybrid Equations (HyEQ) Toolbox in [5D] was used to compute the trajectories. 


4.1 Nominal Case 

The possible solutions to the hybrid system Ha? fall into four categories: always desynchronized, asymp¬ 
totically desynchronized, never desynchronized, and initially synchronized. The following simulation results 
show the evolution of solutions for each category. The parameters used in these simulations are f = 1 and 
e = -0.2. 


4.1.1 Always desynchronized (N G {2,3}) 

A solution toTLN that has initial condition t(0, 0) G A stays desynchronized. Figured shows the evolution of 
such a solution for systems 7^2 and "Hs- Furthermore, as also shown in the figures, for these same solutions, 
the Lyapunov function is initially zero and stays equal to zero as hybrid time goes on. 
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(a) Solutions to H 2 with C 2 = 0.24 and r(0, 0) = (b) Solutions to "Ka with C 2 = 0.32 and t(0, 0) = 
[o,o.i]T £P 2 \X 2 . [0,0.1,0.2]^ eP 3 \X 3 . 




Figure 7: Solutions to'H n that asymptotically converge to the set for TV S { 2 ,3, 7,10}. 


4.1.2 Asymptotically desynchronized {N G {2,3,7,10}) 

A solution of Hat that starts in \ (A U A) asymptotically converges to A, as Theorem 13.41 indicates. 
Figure [7] show solutions to both 7^2 and Ha converging to their respective desynchronization sets. 

For H 2 , if t( 0, 0) = [0,0.1]^, then the initial sublevel set is Ty(c 2 ) with C 2 = 0.24. Using Theorem 13.41 
the time to converge to the sublevel set Lv{ci) with ci = 0.1 leads to M = 7.84. Figure 7(a) shows a 
solution to the system for 10 seconds of flow time. From the figure, it can be seen that V{T{t,j)) « 0.1 at 
{t,j) = (3,4). Then, the property guaranteed by Theorem 13.41 namely, V{T{t,j)) < ci for each (t, j) such 
that t + j > M, is satisfied. Figure [7(b)[ shows a solution and the distance of this solution to A. Notice that 
the initial sub level set is Lv{c 2 ) with C 2 = 0.32. From Theorem 13.41 it follows that the time to converge to 
Ly(ci) with Cl = 0.1 is given by M = 10.14, which is actually already satisfied at (t, j) = (2.2,4). Figure[7] 
show solutions to Hjv that asymptotically desynchronize for N G {7,10}. 


4.1.3 Always Synchronized 

When the impulse-coupled oscillators start from an initial condition r(0, 0) G A, a solution remains in A. 
Figure [ 8 ] shows solutions to H 2 and H 3 that never desynchronize. 

It can be seen that (since T{t,j) G A for all {t,j) G domr) V remains constant. 


4.1.4 Initially Synchronized 

As mentioned in the proof of Theorem 13.31 there exist solutions that are initialized in A and eventually 
become desynchronized. This is due to the set-valuedness of the jump map at such points. Figure [9] shows 
two different solutions to H 2 and H 3 from the same initial conditions r(0, 0) = [0,0, 0]^. Furthermore, notice 
that, for each (t,j), the that Lyapunov function along solutions does not decrease to zero until all states 
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(a) Solutions to H 2 with r(0,0) E A' 2 . (b) Solutions to "Hs with r(0, 0) € ^ 3 . 

Figure 8 : Solutions to Hn that never converge to the set ^ for TV = {2, 3}. 


are non-equal. Recall that from the analysis in Section [3^ when states are equal, the issued solutions are 
outside of the basin of attraction. 


4.2 Perturbed Case 

In this section, we present numerical results to validate the statements in Section [3.41 


4.2.1 Simulations of T-Ln with perturbed jumps 

In this section, we consider a class of perturbations on the jump map and jump set. 

• Perturbation of the threshold in the jump set: We replace the jump set D by Dp := {t :3i £ 
I s.t. Ti = f + Pi} where pi G [0, pi], pi > 0 for each i € I. To avoid maximal solutions that are not complete, 
the flow set C is replaced by Cp := [0, f -b pi] x [0, f -|- P 2 ] x ... x [0, f -b pat]. Furthermore, the components 
of the jump map are also replaced by 


9 Pi {t) = 


{0,Ti(l-be:)} 

(1 + £)Ti 


if Ti = T -b Pi 3j e / \ {*} S.t. Tr = T ■ 
''' < f + Pi 3j € I \ {i} S.t. Tr = f ■ 


if Ti 


( 21 ) 


-Pj 


This case of perturbations is an example of Theorem 13.61 with p affecting only th^jump map. The 
trajectories of the perturbed version of TTat will converge to a region around the set A. Simulations are 
presented in Figures (TUI and ITT] for TV = 2, w = 1, f = 3, and e: = —0.3. 


Figure [TUI shows numerical results for the case when each pi are equal, i.e., pi = P 2 = 0 . 02 . Figure 10 (a) 
shows a solution (solid blue) to the perturbed 7^2 with initial condition t(0,0) = [ 1 . 6 , 2 .!]^ (blue asterisk) 
on the (ti, r 2 )-plane with C (black dashed line), the perturbed flow set Cp (red dashed line), and the 
desynchronization set A (solid green line). From this figure, notice that the solution extends beyond the 
set C and resets at = 3 -b 0.2. The solution converges to a region near the desynchronization set, as 
Theorem 13.61 guarantees. To further clarify the response of 'H 2 to this type of perturbation. Figure |10(b)| 
shows the distance to the set A for 10 solutions with randomly chosen initial conditions t(0, 0) G Cp. Notice 
that for the initial conditions chosen, all solutions converge to a distance of approximately 0.08 by t « 28 
seconds. 


Figure [TT] shows the numerical results for the case when each pi are not equal, i.e., pi yf p 2 . Figure ll(a 
shows 10 solutions from random initial conditions t(0, 0) G Cp with pi = 0.5 and p 2 = 0.4. For this case, the 
solutions converge to a region near A, in that, \T{t,j)\j < 0.22 after approximately 0.28 seconds of flow time. 
Figure [ri(b)| shows 15 solutions when pi = 0.02 and p 2 = 0.01. For this set of simulations, the solutions 
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(a) Solutions to 'H 2 with r(0,0) € X 2 . Notice that 
the solution jumps out of X 2 at (£, j) = (3, 3) and the 
function V begins to decrease after that jump. 




(b) Solutions to Ha with r(0, 0) € A: 3 . At hybrid time 
{t^j) = (1,0) the timer state ri jumps away from the 
other two and begin to desynchronize. At approxi¬ 
mately (tjj') = (4.5,8), all of the states are not equal 
and V begins to decrease. 


Figure 9: Solutions to'H n for N S {2,3} that initially evolve in X and eventually become desynchronized 
due to the set-valuedness of the jump map. 




(a) Solution to H 2 on the (ri, r 2 )-plane with initial (b) Distance to the set A for 10 solutions 
condition r(0,0) = [1.6,2.1]^. with initial conditions randomly chosen from 

[0, X [0,f-|-p2]- The solutions have a dis¬ 

tance that converges to a steady state value of ap¬ 
proximately 0.08 at approximately 28 seconds of 
flow time. 

Figure 10: Solutions to the hybrid system with perturbed threshold, namely, with Up = jr : g 

{1, 2} s.t. Ti=f + Pi} for Pi = p 2 = 0.2. 
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t [seconds] 


(a) Distance to the set A for 10 so¬ 
lutions with random initial conditions 

r(0, 0) € [0, f + pi] X [0, f + P 2 ] with pi = 0.5 
and p 2 = 0.4. 



t [seconds] 


(b) Distance to the set A for 15 so¬ 
lutions with random initial conditions 
t(0, 0) E [0, T + pi] X [0, r + P2] with pi = 0.02 and 
p2 = 0.01. 


Figure 11: Numerical simulations of the perturbed version of I-L 2 with jump set given by Dp = {r : € 

{ 1 , 2 } s.t. Ti = f + Pi} for different values of pi. 


converge to a distance of approximately 0.04 around A after approximately 26 seconds of flow time. These 
simulations validate Theorem 13.61 with p affecting only the jump map, verifying that the smaller the size of 
the perturbation the smaller the steady-state value of the distance to A. 

• Perturbations on the reset component of the jump map: Under the effect of the perturbations 
considered in this case, instead of reseting n to zero, the perturbed jump resets Tt to a value pi € R>o, for 
each i G I. The perturbed hybrid system has the following data: 

/(r) = wl Vr e Cp := C 


and 


Gpir) = bpi(T),...,gpdT)]^ 


Vt G Dp = D 


where, for each i G I, the perturbed jump map is given by 


Pi if Ti = T, Tr<T Vj e / \ {*} 

9i{T) = S [pi^nix + e)} if Ti = f 3j G / \ {i} s.t. Tr = f 
(1 -b e)n if Ti < f 3j G / \ {d s.t. Tr = f 


( 22 ) 


This case of perturbations exemplifies Theorem l3.6l with p affecting only the jump map of Hn- Figures [T^ 
and m show several simulations to this perturbation of Hat. All of the simulations in this section use 
parameters a; = 1, f = 3, e = —0.3, and N = 2. 

The first case of the perturbed jump map Gp considered is for pi = p 2 = 0.02. Figure 12(a) shows a 
solution to the perturbed H 2 from the initial condition t(0, 0) = [2.4, 2.3]^ on the (ti, T 2 )-plane. Notice that 
for T G H such that Ti = f the jump map resets Ti to pi (red dashed line) and not to 0 as in the unperturbed 
case. The solution for this case approaches a region around A, as Theorem 13.61 guarantees. Figure |12(b)| 
shows the distance to the set A over time for 10 solutions of the perturbed system H 2 with initial conditions 
t(0,0) G P 2 \ A 2 . This figure shows that solutions approach a distance of about 0.12 after 25 seconds. 

Now, consider the case where pi d P 2 - Figure [T51 shows the distance to A for two sets of solutions with 
different values for pi and p 2 - More specifically, Figure [l3(a)| shows the case of pi = 0.15 and p 2 = 0.25. 
For this case, it can be seen that the solutions converge after « 28 seconds of flow time and, after that time, 


satisfy |T(t, j)|_^ < 0.25. Figure 13(b) shows the case of pi = 0.02 and p 2 = 0.01. For this case, this figure 
shows that, after ss 28 seconds of flow time, the solutions satisfy |T(t, j)|^ < 0.04. These simulations validate 
Theorem 13.61 with p affecting only the jump map, verifying that the smaller the size of the perturbation the 
smaller the steady-state value of the distance to A. 

• Perturbations on the “bump” component of the jnmp map: In this case, the component 
(1 + e)Ti of the jump map is perturbed, namely, we use t^ = (1 -b e)Ti + Pi{Ti), where pi : R>o —t Pn \ A 
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(a) Solution to "^2 on the (ri, r 2 )-plane with initial 
condition t(0, 0) = [2.4,2.3]”'~. 



initial conditions randomly chosen from C. Most of 
the solutions have a distance that converges to a 
steady state value of approximately 0.12 at about 
25 seconds 


Figure 12: Solutions to the hybrid system %2 with the perturbed jump in (1221) map with pi = p 2 = 0.2. 



(a) Distance to the set A for 10 solutions with ran¬ 
dom initial conditions r(0, 0) € C with pi = 0.15 
and p 2 = 0.25. 



(b) Distance to the set A for 10 solutions with ran¬ 
dom initial conditions r(0, 0) € C with pi = 0.02 
and p 2 = 0.01. 


Figure 13: Solutions to the hybrid system 'H 2 with the perturbed jump map with pi ^ p 2 - 
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(a) Solution "^2 on the (ri, T 2 )-plane with initial (b) Distance to the set ^ for 10 solutions to ^2 with 
condition t(0, 0) = [0.1, 0.2]"''. initial conditions randomly chosen from C. These 

solutions have a distance that converges to a steady 
state value of approximately 0.08 at about 45 sec¬ 
onds. 

Figure 14: Solutions to the hybrid system with perturbed “bump” on the jump map, with pi = P 2 = 0.1. 



(a) Distance to the set A for 10 solutions with ran¬ 
dom initial conditions t(0, 0) G C with pi = 0.15 
and p 2 = 0.1. 



(b) Distance to the set A for 10 solutions with ran¬ 
dom initial conditions t(0, 0) E C with pi = 0.02 
and p 2 = 0.01. 


Figure 15: Numerical simulations of the perturbed version of I-L 2 with the perturbed “bump” on the jump 
map with 'pi ^ p 2 - 


is a continuous function. The perturbed jump map Gp has components Ppi that are given as gi in ([3|) but 
with nil + e) + Pi{Ti) replacing Ti{l+e). 

Consider the case Piin) = piTi with pi G (0, |e|) and let = e + pi G (—1,0). Then reduces to 

= (1 + ei)Ti and the jump map Ppi is given by (HJ with Ei in place of e. This type of perturbation is used 
to verify Theorem 13.61 with p affecting only the “bump” portion of the jump map. Figures fUl and ITS] show 
simulations to 'Hn with the parameters w = 1, r = 3, £ = —0.3, and N = 2. 

Consider the case of H 2 with Gp when pi = p 2 = 0.1, leading to £^i = £2 = 0.2. Figure fUl shows a 
solution on the (ri, T 2 )-plane for this case with initial condition r(0,0) = [0.1, 0.2]^. Notice that the solution 
approaches a region around A (green line), as Theorem 13.61 guarantees. Figure [Td (b) | shows the distance to 
the set A over time for 10 solutions with initial conditions t(0, 0) G C. It shows that solutions approach a 
distance to ^ of « 0.09 after Ri 40 seconds of flow time. 

Next, we consider the case of Gp with £1 ^ 82 - Figure 15(a) shows the distance to A for 10 solutions 
with perturbations given by pi = 0.15 and p 2 = 0.1. For this case, the distance to A satisfies |'r(t, j)|_^ < 0.3 
after Ri 40 seconds of flow time. Figure [T5(b)| shows simulation results with pi = 0.02 and p 2 = 0.01. Notice 
that the smaller the value of the perturbation is, the closer the solutions get to the set A. For this case, 
after Ri 30 seconds of flow time, the distance to A satisfies \T{t,j)\j[ < 0.06. These simulations validate 
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(b) Distance to the set A for 10 solutions of the per¬ 
turbed 1-12 with random initial conditions t(0, 0) € 
C. 


Figure 16: Solutions to the hybrid system 7^2 with perturbed flow map given by the cases covered in 
Section 14.2.21 Figures (a) and (b) show solutions given by the flow perturbation Aw = [0.120,0.134]^ given 
in Section [4.2.21 Note that these figures have a dashed black line denoting the calculated distance from A 
in ([23]) . 


Theorem 13.61 with p affecting only the jump map, verifying that the smaller the size of the perturbation the 
smaller the steady-state value of the distance to A would be. 

4.2.2 Perturbations on the Flow Map 

In this section, we consider a class of perturbations on the flow map. More precisely, consider the case when 
there exists a function {t,j) e-)- c{t,j) such that c{t,j) A c with c as in (fTSl) . Then, from Theorem 13.71 with 
da, we know that 


lim \T{t,j)\j< 


CT 

< 


euj 


euj 


(23) 


Figure [T6| shows a simulation so as to verify this property. The parameters of this simulation are TV = 2, 
w = 1, e = —0.3, f = 4, and Aw = [0.120,0.134]^. It follows from (fT51) that c = 0.0105. Then, from (|17|) . 
it follows that limt_|_j_>oo \T{t,j)\x — 0.1047. Specifically, Figure 16(a) shows a solution on the (ti, T 2 )-plane 
of the perturbed hybrid system 'H 2 with initial condition t( 0,0) = [0,0.01]^. This figure shows the solution 


(blue line) converging to a region around A (between dash-dotted lines about A in green). Figure 16(b) 
shows the distance to the set A of 10 solutions with initial conditions t(0, 0) S C with a dashed line denoting 
the upper bound on the distance in (1231) . Notice that all solutions are within this bound after approximately 
15 seconds of flow time and stay within this region afterwards. 


5 Conclusion 

We have shown that desynchronization in a class of impulse-coupled oscillators is an asymptotically stable 
and robust property. These properties are established within a solid framework for modeling and analysis of 
hybrid systems, which is amenable for the study of synchronization and desynchronization in other impulse- 
coupled oscillators in the literature. The main difficulty in applying these tools lies on the construction 
of a Lyapunov-like quantity certifying asymptotic stability. As we show here, invariance principles can be 
exploited to relax the conditions that those functions have to satisfy, so as to characterize convergence, 
stability, and robustness in the class of systems under study. Future directions of research include the study 
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of nonlinear reset maps, such as those capturing the phase-response curve of spiking neurons, as well as 
impulse-coupled oscillators connected via general graphs. 
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A Appendix 


The following result derives the solution to Tts = 6 with T given in (HOD and b = t1 via Gaussian elimination. 

Lemma A.l For each e G (—1, 0), the solution to Tts = b with T given in (I10[) and b = fl is such that 
its elements, denoted as for each k € {1,2,..., N}, are given by = 


Ef=o (=+i)' 


■T. 


Proof The N x N matrix in (TTOl) and the N x 1 matrix b = t1 leads to the augmented matrix [r|6] given 
by 


■ 1 

0 

0 

0 

0 

0 

r 

0 

(e + 2) 

-(e + 1) 

0 

0 

0 

r 

0 

(e + 1) 

1 

-(e+1) 

0 

0 

f 

0 

(e + 1) 

0 

1 

-(£ + 1) ■■ 

0 

r 

0 

(e + 1) 

0 

0 

0 

■ -(e + 1) 

f 

_ 0 

(e + 1) 

0 

0 

0 

1 

f 


To solve for r^', we apply the Gauss-Jordan elimination technique to (I24p to remove the elements —{e + 1) 
above the diagonal. Starting from the N-th row to remove the — (e + 1) component in the iV — 1 row, and 
continuing up to the second row, gives 


-1 

0 

0 

0 

0 . 

.. 0 

r 

0 


0 

0 

0 . 

.. 0 

YJ^J,\e + iyf 

0 


1 

0 

0 . 

.. 0 

Ef=o'(£+ir^ 

0 


0 

1 

0 

0 

Y^j,\e + iyf 

0 

(e + l)2 + (£ + l) 

0 

0 

0 

0 

r + (1 + e)f 

. 0 

(^+1) 

0 

0 

0 . 

.. 1 

f 


Denoting the augmented matrix in (1251) as [r'|6'], with rf 
element of with k > 2 can be derived from (l24l) as TJ. 2 


^ ( C I 11 ^ 

= T and rf = the solution for each 

rf +Tg = b{ where TJ. 2 denotes the {k, 2) entry 


of r'. Noting that rf can be rewritten as rf 


y.f-Aie+iy^ 

~ —TTT"" 

Ei=o (e+1) 


leads to 


’<r^N-k 

_ 2 ^i=Q 

— v^A/-l 


E iV — 

i=0 



^®For example consider k = 3, the expression reduces to which leads to 

n ^-—1 


E£o^(^ + i)^ 

_ EfEo^(. +1)^ [e[Io'(= +1)^ - (^ +1) EiIo^(. +1)^] ^ + ly 

E[Io'(^ + i)‘ ^ E£o'(s + i)‘^ 
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Lemma A.2 For each x ^ 1, and m, n G N such that n — 1 > m, the finite sum satisfies Y!H=m = 

X — 1 


Proof Let Sn = + 2 ;"*+^ + ... + x” ^ multiply by x to get xSn = 2 ;'"+^ + + 

x"*+3 _|_ _ _|_ 2 ;" Subtracting these expressions leads to 

xSr, -Sn = (a:™+3 + x'"+2 + x'"+3 + ... + a;") - (x’"+3 + x™+2 + x’"+3 + ... + x”-^). 

Then, it follows that Sn = ^ ® 


Lemma A.3 For each x 7 ^ 1, and each m, A G N such that N > m, the finite sum satisfies 


N N-n 


n—m i—0 


xN + (m — A — 2)x + (A — TO + 1) 

(X-1)2 


(26) 


Proof Let Sn = J2n=m ™ 

J2l-im+2) ^ ^;^-(-+3) 2;* + . . . + ^ 

X* leads to 


(1^ . Expanding 5'„ leads to S'™ = + 

i=o 2 ^* + Si=o Then, expanding the sum = x^~'^ + 


A^ —(m+1) N—(m+2) N—(m+3) 1 0 

5„=x'^-'" + 2 ^ F+ ^ x*+ x* + ... + ^x*+^x\ 

2 — 0 2 — 0 2 — 0 2—0 2 — 0 

Expanding leads to 

AT-(m+2) A^—(?72+3) 1 0 

S'„ = 2;^-™ + 2x^-('"+i1+3 x*+ Y a;* + --- + E^* + E®*' 

2—0 2—0 2 — 0 2—0 

The next two sums follow similarly and we arrive to 

N N-n 1 0 

YY.^' = + 4x^-('"+ 3) +... + y] 2 ;* + y] xb 

n—m i—0 i—0 i—0 

Proceeding this way for each sum and noticing that there are exactly (A — to) summations of the form 
Ei=o **> 11 follows then that 

N N—n 0 

+ 2a;^-('"+3) + 3 x^-(™+ 2) + 4x^-(™+3) + ... + (AT _ ,„)2;i + (TV - to + 1) y] x* 

n—m i—0 i—0 

and finally 

Sn = X^-^ + 2x^-1™+3) + 3x^-1™+2) + 4x^-1™+3) + . . . + (A - to)x3 + (A - to + l)x° 


which reduces to 'Yl=m Efco" * m+i^ follows that 

xSn = x^-'"+3 + 2x^-™ + 3 x^-1™+3) + 4x^-1™+2) + ... + (A - to)x2 + (A - to + l)xi 
X^Sn = + 2x^-™+3 + 3x^-™ + dx^-”"-! + . . . + (A - to)x3 + (A - TO + ^X^. 


Then, 

X^Sn - 2xSn + Sn = (x^-™+2 + 2x^-'"+3 + 3x^-™ + dx^-"*"! + . . . + (A - TO 
- 2(x^-”"+3 + 2x^-™ + 3x^-l™+il + 4x^-1™+2) + ... + (A - 
+ (x^-™ + 2x^-1™+3) + 3x^-1™+2) + 4x^-(™+3) + ... + (A ■ 

which reduce to (x—1)^5'„ = x^“™+^+(A—to+1)+(to—A— 2)x, leading to = - 


)x3 + (A — TO + l)x^) 

- to)x^ + (A — to + l)x3) 

— m)x^ + (A — TO + l)x3). 

'^-™+2+(Ar_m+l) + (m-Ar_2)x 

-- 
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